//
// Ephi - simulation of magnetic fields and particles
// Copyright (C) 2007 Indrek Mandre <indrek(at)mare.ee>
// For more information please see http://www.mare.ee/indrek/ephi/
//
// This program is free software; you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation; either version 2 of the License, or
// (at your option) any later version.
//
// This program is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
// GNU General Public License for more details.
//
// You should have received a copy of the GNU General Public License along
// with this program; if not, write to the Free Software Foundation, Inc.,
// 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
//

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <vector>

#include "ephi.hpp"

struct point
{
  vect3d pos;
};

#define EMOD (1e-7*SPEED_OF_LIGHT*SPEED_OF_LIGHT)

class MyE : public StaticElement
{
public:
  MyE (prec_t cpp) :cpp(cpp) { }

  void addFields (const vect3d& pos, vect3d& bfield, vect3d& efield)
  {
    std::vector<prec_t> mags;
    std::vector<size_t> hist;

    for ( size_t i = 0; i < 9; i++ )
      {
        mags.push_back (0.005 + i * 0.005);
        hist.push_back(0);
      }

    for ( size_t i = 0; i < points.size(); i++ )
      {
        prec_t mag = (points[i] - pos).magnitude();
        for ( size_t j = 0; j < mags.size(); j++ )
          {
            if ( mag <= mags[j] * mags[j] )
              {
                hist[j]++;
                break;
              }
          }
      }

    prec_t smooth = 0.05;

    size_t cum = 0;
    for ( size_t j = 0; j < mags.size(); j++ )
      {
        cum += hist[j];
        if ( cum >= 64 )
          {
            smooth = mags[j];
            break;
          }
      }

    vect3d oef, ief;
    oef.clear();
    ief.clear();
    size_t n = 0;
    prec_t tr = 0;
    for ( size_t i = 0; i < points.size(); i++ )
      {
        vect3d r = pos - points[i];
        prec_t rr = r * r;
        prec_t rlen = prec_t_sqrt(rr);
        #define SMOOTH 0.04
        if ( rlen < 1e-23 )
          {
            printf ("Oops!\n");
            return;
          }
        else if ( rlen < smooth )
          {
            tr += rlen;
            ief += (EMOD * cpp / (smooth * smooth * rlen)) * r;
            n++;
          }
        else
          {
            prec_t r3 = rr * rlen;
            oef += (EMOD * cpp / r3) * r;
          }
            
        //if ( ef.magnitude() > 1e12*1e12 )
        //    ef = ef.normal() * 1e12;
      }
    //char b1[128], b2[128];
    //printf ("n=%zu, avg.r=%g, ief=%s, oef=%s\n", n, tr/n, ief.sprint(b1), oef.sprint(b2));
    efield += oef;
    //ief *= 2;
    efield += ief;
  }

  void setCharge (prec_t q)
  {
    cpp = q / points.size();
  }

  void addPotential (const vect3d& pos, prec_t& potential)
  {
  }

  void draw(Screen& screen, const Screen::color& c)
  {
  }

  std::string get_povray_script ()
  {
    return std::string();
  }

  bool inContact (const vect3d& pos)
  {
    return false;
    for ( size_t i = 0; i < points.size(); i++ )
        if ( (pos - points[i]).magnitude() < 0.001 * 0.001 )
            return true;
    return false;
  }

  typedef std::vector<vect3d> points_t;
  points_t points;

  prec_t cpp;
};

#define BFK 1e-7
class MyB : public StaticElement
{
public:
  MyB () { }

  void addFields (const vect3d& pos, vect3d& bfield, vect3d& efield)
  {
    for ( size_t i = 0; i < points.size(); i++ )
      {
        vect3d r = pos - points[i].pos;
        prec_t rr = r * r;
        prec_t rlen = prec_t_sqrt(rr);
        if ( rlen < 1e-23 )
          {
            return;
          }
        else if ( rr < points[i].mlenmlen )
          {
            vect3d bf = (BFK * ELEMENTARY_CHARGE / (points[i].mlenmlen * rlen)) * (points[i].v % r);
            bfield += bf;
          }
        else
          {
            prec_t r3 = rr * rlen;
            vect3d bf = (BFK * ELEMENTARY_CHARGE / r3) * (points[i].v % r);
            //if ( bf.length() > 1.60218e-24 ) bf = bf.normal() * 1.60218e-24;
            bfield += bf;
          }
      }
  }

  void addPotential (const vect3d& pos, prec_t& potential)
  {
  }

  void draw(Screen& screen, const Screen::color& c)
  {
  }

  std::string get_povray_script ()
  {
    return std::string();
  }

  bool inContact (const vect3d& pos)
  {
    //return false;
    for ( size_t i = 0; i < points.size(); i++ )
        if ( (pos - points[i].pos).magnitude() < 0.001 * 0.001 )
            return true;
    return false;
  }

  struct point
  {
    point() { }
    point (const vect3d& pos, const vect3d& v) : pos(pos), v(v) { }
    vect3d pos;
    vect3d v;
    prec_t mlenmlen;
  };

  void recalc ()
  {
    std::vector<prec_t> heap;
    for ( size_t i = 0; i < points.size(); i++ )
      {
        heap.clear();
        for ( size_t j = 0; j < points.size(); j++ )
          {
            if ( j == i ) continue;
            prec_t l = (points[i].pos - points[j].pos).length();
            heap.push_back (l);
            std::push_heap (heap.begin(), heap.end(), std::greater<prec_t>());
            break;
          }
        prec_t m = 0;
        size_t k;
        for ( k = 0; k < 3; k++ )
          {
            std::pop_heap (heap.begin(), heap.end(), std::greater<prec_t>());
            m += heap.back();
            heap.pop_back();
            if ( heap.empty() )
                break;
          }
#if 0
        m /= k;
        m *= 24;
#else
        #define SM 0.05
        m = SM;
        //m *= 3;
#endif
        points[i].mlenmlen = m * m;
      }
    printf ("recalc done!\n");
  }

  typedef std::vector<point> points_t;
  points_t points;
};

class PolywellOctreeBuilder : public OctreeBuilder
{
public:
  PolywellOctreeBuilder (PMap& pmap, const vect3d& c1, const vect3d& c2, size_t limitMemoryMB = 100)
    : OctreeBuilder (c1, c2, 1, 2, limitMemoryMB, 0.05), pmap(pmap)
  {
  }

  class mydifferator : public differator<7>
  {
    public:
      mydifferator(PMap& pmap) : differator<7>(vect3d(1e-04, 1e-04, 1e-04)), pmap(pmap) { }

      void getValues (const vect3d& pos, prec_t *vals)
      {
        CubicScalar pot;
        CubicVector cef, cbf;
        pot.clear();
        cef.clear();
        cbf.clear();
        pmap.addCubicValues (pos, pot, cbf, cef);
        vals[4] = cef.f.x;
        vals[5] = cef.f.y;
        vals[6] = cef.f.z;
      }

    private:
      PMap& pmap;
  };

  void addValues (const vect3d& pos, CubicScalar *scalars, CubicVector *vectors)
  {
#if 1
    pmap.addCubicValues (pos, scalars[0], vectors[0], vectors[1]);
#else
    memset (dfdx, 0, 7 * sizeof (prec_t));
    memset (dfdy, 0, 7 * sizeof (prec_t));
    memset (dfdz, 0, 7 * sizeof (prec_t));
    memset (d2fdxdy, 0, 7 * sizeof (prec_t));
    memset (d2fdxdz, 0, 7 * sizeof (prec_t));
    memset (d2fdydz, 0, 7 * sizeof (prec_t));
    memset (d3fdxdydz, 0, 7 * sizeof (prec_t));

    mydifferator diff(pmap);
    diff.getValues (pos, f);
    diff.getDifference (pos, dfdx, dfdy, dfdz, d2fdxdy, d2fdxdz, d2fdydz, d3fdxdydz);
#endif
  }

  bool inContact (const vect3d& pos)
  {
    return false;
  }

private:
  PMap& pmap;
};

int main (int argc, char *argv[])
{
{
  PMap pmap(vect3d(0.5, 0.5, 0.5), 0.0075, 0, 0.1);
  pmap.addQ (vect3d(0.00001, 0.00001, 0.00001), vect3d(1, 0, 0), 1);
  pmap.recalc();
  Octree tree;
  PolywellOctreeBuilder builder(pmap, vect3d(-0.5, -0.5, -0.5), vect3d(0.5, 0.5, 0.5));
  builder.build (tree);
  Statics statics;
  statics.addStaticElement (new StaticCubePolywellOctree (tree, true, true, false, 0.15, 0, 0.01, 1, 1));
  ElectroDynamics dyn(statics);
  Scene scene(dyn);
  Screen screen(800, 800, 0.8, 0.8);
  screen.set_axis (vect3d(0, 0, 0), vect3d(1, 1, 0).normal(), vect3d(0, 0, 1));
  scene.calc (screen, Scene::CALC_EFIELD);
  scene.set_coloring (0, true);
  scene.render_map (screen);
  screen.write ("tt2.bmp");
  scene.render_lic (screen);
  screen.write ("tt2_lic.bmp");
  exit(0);
}

  prec_t fmin, fmax;
  if(1)
  {
  Statics statics;
  ElectroDynamics dyn(statics);
  statics.addStaticElement (new StaticUniformSphere (vect3d(0, 0, 0), 0.1, 1.0));
  Screen screen(800, 800, 0.75, 0.75);
  Scene scene(dyn);
  scene.calc (screen, Scene::CALC_EFIELD);
  fmin = scene.get_fmin();
  fmax = scene.get_fmax();
  scene.set_coloring (32, true);
  scene.render_map (screen);
  screen.write ("tt2_us_o.bmp");
  scene.set_coloring (0, true);
  scene.render_lic (screen);
  screen.write ("tt2_us_lic_o.bmp");
  }

  if (1)
  {
  Statics statics;
  ElectroDynamics dyn(statics);

  MyE *myq = new MyE (0);

  #define RADIUS 0.10

  for ( size_t i = 0; i < 30000; i++ )
    {
      prec_t x = prec_t_drand() * (2 * RADIUS) - RADIUS;
      prec_t y = prec_t_drand() * (2 * RADIUS) - RADIUS;
      prec_t z = prec_t_drand() * (2 * RADIUS) - RADIUS;
      vect3d pos(x, y, z);
      if ( pos.magnitude() > 0.1 * 0.1 )
          continue;
      myq->points.push_back (pos);
    }

  myq->setCharge (1);
  statics.addStaticElement (myq);

  Statics s2;
  s2.addStaticElement (new StaticUniformSphere (vect3d(0, 0, 0), 0.1, 1.0));

  for ( size_t i = 0; i < 20; i++ )
    {
      vect3d ef1, ef2, tmp;
      vect3d pos(-(prec_t)i * 0.01, 0, 0);
      statics.getFields (pos, tmp, ef1);
      s2.getFields (pos, tmp, ef2);
      char b0[128], b1[128], b2[128];
      size_t c = 0;
      for ( size_t j = 0; j < myq->points.size(); j++ )
        {
          if ( (myq->points[j] - pos).magnitude() < SMOOTH * SMOOTH )
              c++;
        }
      prec_t err = (ef1-ef2).length() / ef1.length();
      printf ("%s: simu=%s, real=%s error=%.2f%%(close=%u)\n", pos.sprint(b0), ef1.sprint(b1), ef2.sprint(b2), 100 * err, (unsigned int)c);
    }
  //exit(0);

  Screen screen(800, 800, 0.75, 0.75);
  Scene scene(dyn);
  scene.calc (screen, Scene::CALC_EFIELD);
  scene.calibrate_field (fmin, fmax);
  scene.set_coloring (32, true);
  scene.render_map (screen);
  screen.write ("tt2_us_myq.bmp");
  scene.set_coloring (0, true);
  scene.render_lic (screen);
  screen.write ("tt2_us_lic_myq.bmp");
  }

  if (0)
  {
  Statics statics;
  ElectroDynamics dyn(statics);

  MyB *myb = new MyB();

#if 0
  for ( size_t i = 0; i <= 40; i++ )
    {
      for ( size_t j = 0; j <= 40; j++ )
        {
            {
              prec_t x = ((prec_t)i - 20) * 0.05;
              prec_t y = ((prec_t)j - 20) * 0.05;
              vect3d pos(x, y, 0);
              if ( pos.length() > 0.31 )
                  continue;
              //pos.print();
              myb->points.push_back (MyB::point(pos, vect3d(0, 0, 1)));
            }
        }
    }
#endif
  for ( size_t i = 0; i < 5000; i++ )
    {
      prec_t x = (prec_t_drand() - 0.5) * 1;
      prec_t y = (prec_t_drand() - 0.5) * 1;
      vect3d pos(x, y, 0);
      if ( pos.length() > 0.31 )
        {
          i--;
          continue;
        }
      //pos.print();
      myb->points.push_back (MyB::point(pos, vect3d(0, 0, 1)));
    }
  myb->recalc();
  statics.addStaticElement (myb);

  Statics s2;
  s2.addStaticElement (new StaticUniformSphere (vect3d(0, 0, 0), 0.31, 1.0));

  Screen screen(800, 800, 0.75, 0.75);
  //screen.set_axis (Screen::PLANE_XY, 0.01);
  Scene scene(dyn);
  scene.calc (screen, Scene::CALC_BFIELD);
  scene.set_coloring (32, true);
  scene.render_map (screen);
  screen.draw_cross (vect3d(-0.31, 0, 0), Screen::BLUE);
  screen.draw_cross (vect3d(-0.31-SM, 0, 0), Screen::BLUE);
  screen.write ("tt2_cmc_myb.bmp");
  scene.set_coloring (0, true);
  scene.render_lic (screen);
  screen.write ("tt2_cmc_lic_myb.bmp");
  }

  return 0;
}

